R Notebook

library(tidyverse)
library(taxonomizr)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

genome_data <- read_csv("../data/mmetsp_ncbi_genome_info.csv") %>%
  mutate_at("sample_id", as.character)

── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  .default = col_character(),
  sample_id = col_double(),
  `GenBank BioSample` = col_double(),
  project_id = col_double(),
  latitude = col_double(),
  longitude = col_double(),
  taxon_id = col_double(),
  depth = col_double(),
  day_portion_of_day_night_cycle_in_hours = col_double(),
  light = col_double(),
  experimental_salinity = col_double(),
  experimental_temperature = col_double(),
  night_portion_of_day_night_cycle_in_hours = col_double(),
  ph = col_double(),
  environmental_salinity = col_double(),
  collection_time = col_time(format = ""),
  environmental_temperature = col_double(),
  phosphate = col_double(),
  nitrate = col_double(),
  iron = col_double(),
  trace_elements = col_double()
  # ... with 21 more columns
)
ℹ Use `spec()` for the full column specifications.

genome_data

report_filenames <- list.files("../results/krakenReports/")

report <- NULL

for(filename in report_filenames){
  data <- read_tsv(file.path("../results/krakenReports", filename),
                   col_names = c("percent", "clade_count", "taxon_count", "kmers", "unique_kmers", "level", "taxid", "name" ))
  
  data$file <- filename
  
  
  report <- bind_rows(report, data)
  
}

── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)


── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  percent = col_double(),
  clade_count = col_double(),
  taxon_count = col_double(),
  kmers = col_double(),
  unique_kmers = col_double(),
  level = col_character(),
  taxid = col_double(),
  name = col_character()
)

report <- report %>%
  mutate(duplicity = kmers/unique_kmers,
         ratio_unique = unique_kmers/kmers,
         taxid = as.character(taxid)) %>%
  separate(file, c("sample_id", "sample_name_main", "fullSpeciesName"), sep = "_-_", remove = F) %>%
  filter(clade_count > 1)

taxonomyInfo <- getTaxonomy(unique(report$taxid), sqlFile = "../data/accessionTaxa.sql") %>%
  data.frame(taxid = row.names(.), row.names = NULL) %>%
  # remove spaces in taxid
  mutate(taxid = str_trim(taxid))

report <- report %>%
  left_join(., taxonomyInfo, by = "taxid") %>%
  left_join(., genome_data[, c("sample_id", "taxon_id")], by = "sample_id") %>%
  left_join(., genome_data[, c("sample_id", "species_taxid")], by = "sample_id") %>%
  rename(target_taxid = taxon_id,
         target_species_taxid = species_taxid) %>%
  mutate(is_target_species = (taxid == target_taxid | taxid == target_species_taxid))

report

report %>%
  filter(superkingdom == "Bacteria")

nrow(report)

[1] 70584

sum(report$is_target_species)

[1] 24

taxonomyInfo <- getTaxonomy(unique(report$taxid), sqlFile = "../data/accessionTaxa.sql") %>%
  data.frame(taxid = row.names(.), row.names = NULL) %>%
  # remove spaces in taxid
  mutate(taxid = str_trim(taxid))

taxonomyInfo

write_csv(report, "../data/duplicity_report.csv")

report %>%
  filter(superkingdom == "Bacteria")

report %>%
  filter(superkingdom == "Bacteria",
         grepl("^(S|G).*", level)) %>%
  ggplot(aes(x = taxon_count/unique_kmers))+
  geom_histogram(bins = 30)+
  #geom_density()+
  xlim(0, 10)

report %>%
  filter(
         is_target_species) %>%
  ggplot(aes(x = log(taxon_count)))+
  geom_histogram(bins = 30)#+

  #geom_density()+
  #xlim(0, 10)

report %>%
  filter(superkingdom == "Bacteria",
         grepl("^(S|G).*", level),
         taxon_count != 0) %>%
  ggplot(aes(x =log10(taxon_count)))+
  geom_histogram(bins = 30)+
  #geom_density()+
  xlim(0, 10)+
  ylim(0,5000)

report %>%
  filter(superkingdom == "Bacteria",
         taxon_count/unique_kmers != 0) %>%
  ggplot(aes(x = log(taxon_count/unique_kmers)))+
  geom_histogram(bins = 30)+
  #geom_density()+
  xlim(0, 10)

report %>%
  filter(superkingdom == "Bacteria",
         taxon_count/unique_kmers != 0) %>%
  ggplot(aes(x = (taxon_count/unique_kmers)))+
  geom_histogram(bins = 30)+
  #geom_density()+
  xlim(0, 100)+
  ylim(0, 5000)

report %>%
  filter(superkingdom == "Bacteria") %>%
  ggplot(aes(x = log(clade_count)))+
  geom_histogram(bins = 30)+
  #geom_density()+
  xlim(0, 10)

report %>%
  filter(!is_target_species,
         duplicity > 0) %>%
  ggplot(aes(x = (duplicity)))+
  geom_histogram()+
  scale_x_log10()

report %>%
  #filter(is_target_species) %>%
  ggplot(aes(x = duplicity))+
  facet_wrap(~ifelse(is_target_species,"Target Species", "Other Species"), scales = "free_y")+
  geom_histogram()+
  scale_x_continuous(labels = scales::number_format(accuracy = 1),
                     trans = "log10",
                     limits = c(1, 1000),
                     n.breaks = 8)+
  theme_bw()

report %>%
  #filter(is_target_species) %>%
  ggplot(aes(x = taxon_count))+
  facet_wrap(~ifelse(is_target_species,"Target Species", "Other Species"), scales = "free_y", "free_x")+
  geom_histogram()+
  scale_x_continuous(labels = scales::number_format(accuracy = 1),
                     trans = "log10")+
  theme_bw()

Coercing `nrow` to be an integer.NAs introduced by coercion`nrow` is missing or less than 1 and will be treated as NULL.

report %>%
  filter(!is_target_species) %>%
  ggplot(aes(x = taxon_count))+
  geom_histogram()+
    scale_x_continuous(labels = scales::number_format(accuracy = 1),
                     #trans = "log10",
                     n.breaks = 10,
                     limits = c(1, 500))+
  theme_bw()

report %>%
  filter(superkingdom == "Bacteria") %>%
  ggplot(aes(x = (duplicity)))+
  geom_histogram(bins = 45)+
  #geom_density()+
  xlim(0, 50)

report %>%
  filter(superkingdom == "Bacteria") %>%
  ggplot(aes(x = log(unique_kmers)))+
  geom_histogram(bins = 30) +
  xlim(0, 10)

exp(2.5)

report %>%
  filter(superkingdom == "Bacteria") %>%
  ggplot(aes(x = log(unique_kmers)))+
  geom_histogram(bins = 30) +
  xlim(0, 10)

report %>%
  ggplot(aes(x = log(unique_kmers)))+
  geom_histogram(bins = 50)# +
  xlim(0,25)

nrow(report)

unique(report$level)

report %>%
  filter(clade_count > 200,
         taxon_count > 10,
         unique_kmers > 10,
         duplicity < 70
         )

report %>%
  filter(superkingdom %in% c("Archaea", "Bacteria"),
         clade_count > 1000,
         taxon_count > 500,
         unique_kmers > 10,
         duplicity < 10,
         grepl("^(S|G).*", level)
         #level %in% c("G","G1", "G2", "S", "S1", "S2", "S3")
         )

report %>%
  ggplot(aes(x = duplicity))+
  geom_histogram(bins = 30) +
  #geom_density()+
  xlim(0, 100)

report %>%
  ggplot(aes(x = ratio_unique))+
  geom_histogram(bins = 30)